AN EXACT SOLUTION TO THE TEMPERATURE EQUATION 
IN A COLUMN OF ICE AND BEDROCK 



ED BUELER 



1. The problem 

The goal here is fairly straightforward. We want a solution of a pure conduc- 
tion problem in ice and bedrock. This solution needs to be suitable for verifying a 
numerical scheme for conservation of energy. This solution will also help with the 
construction of an approximate polythermal scheme. We will use this exact solu- 
tion in the context of a coupled ice flow and conservation of energy model, namely 
PISM [9]. This exact solution will form one of a suite of verification tests for PISM 

In particular, we will find a function T{z, t) with the following properties 

TiH,t)=T,, 

dT d'^T 
Pici-g^ = ki-g^ {0<z< H), 

T(0+,t)=T(0-,t), 
dT dT 

dT d^T 
PrCr-q^ = kR-^ {-B<z< 0), 

dT 

-kR—{-B,t) = G. 
dz 

The two conditions at the ice/rock interface z = are continuity of temperature and 
of heat flux, respectively. 

The ice thickness is H > and the bed thickness is 5 > 0; representative values 
used here are 

5 = 1000 m and iJ = 3000 m. 

The ice occupies < z < H and has density p/, specific heat capacity c/, and 
conductivity kj. The bedrock occupies —B < z < and has density pr, specific heat 
capacity cr, and conductivity k^. Reasonable values of these constants are given in 
the C implementation at the end. The constant value Tg of the surface temperature 
will be 223.15 K or —50 °C. The value of the geothermal flux used here is 

G = 42 mW/m^. 
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Let US take as our initial condition an (absolute) temperature which is a linearly- 
increasing function of the depth below the surface of the ice: 

(1) T{z,0) = Ts + (l){H - z), T, = 223.15 K, = 0.0125 K m"^ 

Figure [1] includes a graph of this linear initial condition, which warms from —50 °C 
at the surface to °C at the base of the bedrock layer (i.e. at depth 1000 m into the 
bedrock). 



30 




z 



Figure 1. Final temperature T(z, +oo) (solid) and initial tempera- 
ture T{z,0) (dashed). Within the ice the temperature must actually 
remain below the pressure- melting temperature Tp.^p{z) (dash-dotted). 

In fact we will slightly revise this initial condition in Section [31 In particular, for 
numerical accuracy reasons, it will be desirable to use an initial condition with a finite 
eigenfunction expansion. The graph of the initial condition in Figure [1] is accurate 
at printer resolution, however. 

As noted, a goal is to verify parts of the thermomechanical model in PISM. On 
the other hand, PISM is primarily a three (spatial) dimensional model for the flow 
of ice, coupled with the thermodynamics of the ice and the bedrock. Therefore, in 
using the exact temperature solution here for verification, we will suppose that the 
conditions for the full, coupled model are ice of constant thickness H everywhere, 
accumulation which is identically zero, and a fiat bed. Then PISM will predict no 
flow. In particular, the advection, strain-heating, and basal frictional heating parts of 
the general conservation of energy equations are each identically zero. So in this case 
we see that the temperature problem above is all that remains to solve in the full, 
coupled model. Note that other tests fully verify the conservation of energy numerical 
scheme in flowing ice [21 [3] , but they do not include heat storage in bedrock. 
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We have not included melting in the above. Recall that the pressure-melting 
temperature in the ice is Tp^p{z) = Tq — l3{H — z). (We suppose Tq = 273.15 K and 
13 = 8.66 X 10-^K/m for concreteness. With these constants Tpmp(O) = -2.598 °C.) 
For sufficiently large t > 0, the solution to the above problem has T{z,t) > Tpinp(-2) 
for some locations z > 0. At such locations the above model no longer applies because 
there will be partial melting in the ice. For verification purposes we are interested in 
the ffist time at which melting occurs. 

2. Finding an eigenfunction expansion 

We find a classical kind of solution to this classical kind of problem. First we 
transform our inhomogeneous problem to a homogeneous one. Let 

(z/kj-H/kr, 0<z<H 
[z/kn - H/ki, -B <z<0. 

Define the rescaled temperature 

e{z,t) =T{z,t)-Ts + GP{z). 

It is straightforward to check that T{z,t) solves the original problem if and only if 
9{z, t) solves 

9{H,t) = 0, 
PicA = kiO,, {0<z <H), 

M.(o+,t) = M.(o-,t), 

pRCR0t = kRO,, {-B <z <0), 

e,{-B,t) = o. 

This boundary value problem is linear and homogeneous. (Note we have also switched 
to subscript notation for derivatives.) 

The rescaled temperature 6 has initial condition 

(3) e{z,0) = GP{z) + (j){H-z). 

We expect that the above problem for 9{z,t) is well-posed \7]. Furthermore we 
expect that limf^+oo 0{z, t) = 0, that is, we expect that the problem is asymptotically 
stable. Thus we expect 

T{z,+oo) = Ts-GP{z). 

Of course this would violate the requirement that T < Tpmp within the ice. Nonethe- 
less this final state is worth graphing along with the initial state and the pressure- 
melting temperature, as in Figure [H 
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Next we separate variables and seek eigenfunctions. Preliminary thoughts might 
go like this: If 9{z,t) = fj{z)g{t) on the interval < z < H and 6{z,t) = fn{z)g{t) 
on the interval —B < z < then we have 

PjCj- = ki— and PrCr- = kR—. 
9 Ji 9 jR 

This separated form for 6{z,t) must have, and does have, the same dependence on t 
in both the ice and the bedrock. The solution must satisfy boundary conditions at 
z = H and z = —B. As usual for the heat equation, the solution decays exponentially 
in time and is (roughly) sinusoidal in space. The conditions at z = correspond to 
continuity of the solution and of the heat flux, and this means a change in amplitude 
for the sinusoid because the conductivity changes. 

A conclusion to the above thoughts is an ansatz for separated solutions^ 



(4) 9{z,t) = e 



sm{a{H - z)), < z < H, 
jcos{(3{B + z)), -B<z<0. 



The eigenvalues are denoted A. We will see that they form a countable sequence 
< Ao < Ai < A2 < . . . which tends to positive infinity. The eigenfunctions are the 
spatial parts of the corresponding ansatz solutions. 

The constants A, a, f3, 7 are determined by the connection conditions and the PDEs 
themselves: 

(5) P/C/A = kja'^, 

(6) sm{aH) =-fcos{pB), 

(7) akj cos{aH) = P^ykRsin^PB), 

(8) PrCrX = kR(3\ 
Conditions ([5]) and ([8]) combine to eliminate A and give 

(9) l3 = Za 

where 

' prCr ki 
kR pici 

On the other hand, conditions (EI) and (JTj) combine to eliminate 7. Indeed, using 
as well, after clearing fractions one gets 

(10) Asm{Ha.) sm{ZBa.) = cos{Ha) cos{ZBa) 
where 

A = ^Z. 

ki 



^The reader who does not like this language may confirm that, at the end, we have a full spectral 
resolution of our discrete spectrum, self-adjoint operator. 
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Using trigonometric identities one can rewrite (jlOl) as 



'111 



A 



cos{{H - ZB)a) = cos{{H + ZB)a) 



Note that < A < 1 for reasonable values of density, specific heat capacity, and 
conductivity for ice and bedrock. Thus, 



A 



A + l 



< 1, 



so fill I) equates two sinusoidal functions, with the left-hand function of smaller mag- 
nitude and lower frequency. 

We have arrived at a visualizable stage. Equation (fTTl) determines countably many 
discrete values a = > 0, A; = 0, 1, 2, . . . , as shown in Figure [2l Solutions of (fTTl) 
must occur between each consecutive extrema of the higher amplitude and higher 
frequency cosine on the right side of the equation. Though (fTTl) is transcendental, 
accurate solutions are easily found by good numerical methods like Brent's method 
[6]. In particular, one can bracket each solution, and Brent's method maintains such 
a bracket as it converges to a root. 

We need only find positive solutions a of equation (fTTj) . They will form a positive 
increasing sequence < ao < «i < «2 < • • • • Also, as a special case which may 
be used to check formulas, note that if the material constants p, c, and k are non- 
physically assumed to be the same for ice and for bedrock then Z = 1 and A = 1 so the 
equation we solve is just cos((i/+ 5) a) = 0. In this case = ((2A; + l)7r)/(2(if+i?)). 




0.006 



Figure 2. A picture of equation (fTTl) . There is exactly one solution 
ttfe > per half-cycle of the higher amplitude cosine. 



6 



ED BUELER 



Once we ak then from ([9]) we get a corresponding sequence [3k. From (EI) or (JTj) 
we find 7^. From ([5]) or ([H]) we get the (positive) eigenvalues themselves. In fact, 
using the constants specified below in the C implementation, we get the spectrum 
{Afc} shown in Figure [3l 
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Figure 3. The first 30 eigenvalues Afc, /c = 0, 1, . . . , 29. 



The unnormalized eigenfunctions are 



9k.{z) 



sin(afc(if - 2;)), < 2; < if , 
7fcCos(/?fe(5 + 2;)), -B<z<Q. 



Those 9k corresponding to the five smallest (most important) eigenvalues A^ are 
shown in Figure HI 

The eigenfunctions 9k{z) are solutions of a Sturm-Liouville problem [1]. Thus 
they are an orthogonal set with respect to an appropriate inner product. This inner 
product includes the coefficients used in computing the thermal energy. In fact, recall 
that jjjypcT dxdydz is the internal (specific) heat energy stored in a material with 
temperature T occupying a volume V. So, if f{z),g{z) are integrable functions on 
—B<z<H,we define the inner product: 



:i2) 



"0 i-H 

{f,9)-=pRCRl fiz)g{z)dz + pici f{z)g{z)dz. 

-B Jo 
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Figure 4. Unnormalized eigenfunctions 6o{z), . . . ,94{z). Note a 
change in amplitude at 2; = 0. 



An easy calculation computes inner products of the (as yet) unnormalized eigen- 
functions, as follows. First we transform to doable integrals, 



^ pRCRlkli [ cos{Pk{B + z)) cos{Pi{B + z)) dz 

J-B 

+ PiCi / sin{ak{H — z)) sm{ai{H — z)) dz 
Jo 

nB nH 

= PRCRlkli / cos{(3kx) cos{Pix) dx + picj I sin(Q;fey) sin{aiy) dy. 
Jo Jo 

Now there are two cases, li k — I then we have a formula for normalization constants: 



H 



Ok, Ok ) = ]^PrCriI 1 + cos(2/3fcx) dx + ^pid / 1 - cos(2afcj/) dy 



- 2^-'-^^ + J + 2^^^^ 2^ 

= J {prCr^IB + piCiH) + PrCr^I siniPkB) cos(/3feS) 



2ak 



PiCj sm{akH) cos{akH). 



ED BUELER 



This expression simplifies further using the properties of the eigenf unctions: 

^k = 7; {prCriIB + piCjH) + ^ pRCR-flsm{f3kB) cos{(3kB) 
Z Zpk 

pici-ik cos{fJkB) ; — sm{(jkB) 



2a k akkj 

7^ {prCriIB + piCiH) + ^„ sin(/3fc5) cos{(3kB) (pRCRalkj - pidPlkR) 
I zpkaf^Ki 

^ {prCriIB + piCiH) . 

The starred equahty follows from equations IQ and ([7]). The double-starred equality 
follows from equations ([5]) and ([8]). 

If k I we get (^9k, 9iJ = 0, but we omit the details. 

Thus the normalized eigenf unctions are 

^,(^) = = 1. U'^^iMH - z)), 0<z<H, 

Xk Xk 1^7^ cos{Pk{B + z)), -B < z < 0. 

3. The solution to the time-dependent problem 
The solution to the time-dependent problem for 0{z,t) is the infinite series 

oo 

(13) 9{z,t) = Y,Cke'^'''9k{z) 

k=0 

where 9k{z) are the normalized eigenfunctions computed above, and Ck = {Ok, 9{t = 0)). 
In fact, 

Ck = pici / 9kiz) {GP{z) + (t){H - z)) dz + pRCR / 9kiz) {GP{z) + 0(/7 - z)) dz 

Jo J-B 

from equation ([3]). Also note P{z) is given in equation ([2]). We can naturally describe 
Ck as a linear combination of definite integrals: 

(14) Ck = {pjci II + pRCR^k IT) , 
where 

II = sm{ak{H - z)) (c (^^ - £^ + 0(if - z)^ dz 



and 



-'fc 



y"°^ cos{Pk{B + z)) (c - ^) + - '^)) d^- 



These are elementary integrals, though it is easy to get things wrong anyway. They 
simplify to 

(15) Ik = - ^ [sin(afeiJ) - (akH) cos(afcif)] , 
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(16) II = - 0^ [cos{hB) - 1 + {hB) sm(/3feS)] 

The temperature itself (not rescaled) is given by 

(17) T{z,t) = e{z,t)+Ts-GP{z). 

Formulas f|T3|) . (|T5|) . f|T6l) . and f|T71) together form the time-dependent solution 
to the initial value problem specified so far. 

Now, the infinite sum converges quickly for large times but it converges rather 
slowly for t = 0. This fact relates to the poor differentiability of the initial state 
(times the diffusivity, that is), and it is a common situation for conduction problems 
[1]. To avoid any concern with convergence at t = 0, we redefine the initial state 
to have a finite eigenfunction expansion. That is, we replace equation ([3]) with the 
revised condition 

29 

(18) e{z,o) = Y,Ckek{z) 

k=0 

where the coefficients Ck are given exactly as before by equations dH]), ( 0^51) . and ( fT6l) . 
This represents a change of the initial condition by a maximum of only about 0.001 
K, so to printer or screen accuracy this is not important, and indeed the upper limit 
of the sum = 29 was chosen for such reasons. But that detail is not important. 
Rather, the point is that by making this change any concerns about evaluating the 
exact solution to high accuracy are immediately resolved, and this is our goal. Note 
that this change also means that the time-dependent solution has a finite expansion: 

29 

(19) T{z, t)=Ts- GP{z) + J2 Cke'''^%{z). 

k=0 

4. Verification of PISM using this exact solution 

The exact solution given by equation (fT9|) is verification Test K in PISM [9]. As 
previously noted, PISM is a three-dimensional ice fiow simulation program which 
includes many coupled physical models. Here we use Test K to verify the part of 
PISM which relates to the simulation of heat conduction. That is, PISM contains a 
semi-implicit finite difference approximation of a shallow (continuum) approximation 
of the conservation of energy equation. Our use of Test K for verification concerns 
only the pure conduction aspect of that scheme. 

We note that bugs can and have appeared in the part of PISM which numerically 
approximates the point in the bedrock where the geothermal fiux is applied and 
at the switch of material properties from ice to bedrock. Of course in a many- 
physical-models code like PISM there are many contributions to the approximation 
of conservation of energy at the ice-bedrock interface, including basal melting and 
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frictional heating, and thus the numerical scheme for grid points at the base of the 
ice is comphcated. An exact solution is helpful for debugging such details even if it 
only verifies a sub-model of the full "multi-physics" . 

The vertical grid in PISM has, for now, constant spacing Az, and indeed this 
spacing is equal in both the ice and the bedrockl^ 

The numerical scheme in PISM for the energy equation is documented in the Ap- 
pendices of [Sj. The scheme is semi-imphcit generally, but when restricted to pure 
conduction in a column of ice, as here, it is fully-implicit. That is, it corresponds to 
centered-spatial-differencing and backward Euler method in time and thus it has local 
truncation error 0(At, A^;^). It is unconditionally stable (for pure conduction). In- 
deed, the numerical issues associated to advection and to strain heating, as discussed 
in [3], are not important here. 

To verify using Test K we choose a refinement path [8] with /\z = 100, 50, 25, 12.5, 
and 6.25 meters. As long as At is reduced appropriately, which means At = CAz^ 
for some appropriate C, this gives a refinement path along which the error should 
decay by a factor of four at each refinement. In fact we use At = 400, 100, 25, 6.25, 
and 1.5625 years, so in fact C = 0.04. Because the exact and numerical solutions 
have constant dependence on x and y, the horizontal grid is fixed as at a convenient 
(very coarse) level. 

As shown in Figures [5] and [6], which are admittedly boring figures, the maximum 
and average numerical errors at all points within the ice and within the bedrock do 
decay to zeroj^ As shown in the figures, fitting the average error versus Az to a curve 
of the form (err) = A{AzY gives r = 2.01 for the approximation within the ice and 
r = 2.00 within the bedrock. 

This suggests that the numerical scheme is achieving the optimal rate, that is, the 
local truncation error is reflected in the global approximation error. 

Note that along this refinement path, as Az is reduced by a factor of two we must 
reduce At by a factor of four if we want the time part of the local truncation error to 
contribute a comparable fraction of the error. Along this refinement path the amount 
of computational work per step therefore goes up by a factor of two but the amount of 
computational work per model year goes up by a factor of eight. This statement turns 
out to be slightly pessimistic, because Figure [7] suggests that, running in parallel with 
two processors, the run time for PISM is related to the —2.5 power of Az. That is, 
instead of a halving of Az generating a slowdown of a factor of 8 = 2^, there seems 
to be a slowdown by a factor of only 2^-^. This is probably related to the increasing 
efficiency of the code as more points are computed in each column. 



^This statement applies to PISM in October 2007, but future versions may be change. Such 
changes to the grid are exactly the kind of numerical issue which motivates building and documenting 
a suite of exact solutions for verification. 

■^Mathematical readers should note that we are reporting both L°° and error. 
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Figure 5. Maximum (squares) and average (circles) errors made by 
PISM in approximating the temperature within the ice in Test K. 
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Figure 6. Maximum (squares) and average (circles) errors made by 
PISM in approximating the temperature within the bedrock in Test K. 

Finally, the critical time t when ^(O, t) first exceeds pressure-melting is between 
133, 000 years and 134, 000 years. Indeed, by bisection on the exact solution, it 
must be within a year of 133,465 years. With a modestly refined grid with Az = 
25 m we see the numerical approximation first has T(0, 1) reach pressure melting 
between 133, 470 and 133, 480 model years. This seems close enough, and no further 
verification has been pursued. 
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Figure 7. Run time for PISM to complete Test K using two processors. 
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Appendix A. Reference implementation of Test K 

This Appendix contains a C code which accepts t and z and computes the (abso- 
lute) temperature T given by equation ( |T9i) . That is, this code evaluates Test K in 
PISM. It has only been compiled with the GNU gcc compiler, and the reader may 
note that it is not particularly written for efficiency or speed. 

The file which contains the code is called exactTestK. c, and it is listed verbatim. 
A header file exactTestK. h exists in the PISM source tree, but listing it here would 
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add no information so it is omitted. Likewise there is also a simple example program 
simpleK. c for evaluating the exact solution, but we do not list it. 

The procedure exactKO in exactTestK.c is devoted to evaluating the exact so- 
lution using saved values of ak- These values may be recomputed using the part 
of the code which is delimited by "#if COMPUTE. ALPHA" and "#endif " . This latter 
part uses Brent's method, as implemented in the (GNU Scientific Library , to solve 
equation (fTT]) numerically to about 14 digits of accuracy (in double precision). 

The numerical approximation of conservation of energy within PISM is, of course, 
not listed here. The latest revision can be found at the PISM download site 191. 



Copyright (C) 2007 Ed Busier 
This file is part of PISM. 

PISM is free software; you can redistribute it and/or modify it under the 
terms of the GNU General Public License as published by the Free Software 
Foundation; either version 2 of the License, or (at your option) any later 
version. 

PISM is distributed in the hope that it will be useful, but WITHOUT ANY 
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 
details . 

You should have received a copy of the GNU General Public License 
along with PISM; if not, write to the Free Software 

Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 

*/ 

#include <stdio.h> 
#include <math.h> 
#include <gsl/gsl_errno.h> 
#include <gsl/gsl_math.h> 
#include <gsl/gsl_roots .h> 
#include "exactTestK.h" 



#def ine 


pi 


3.1415926535897931 




#def ine 


SperA 


31556926.0 


/* 


seconds per year; 365.2422 days */ 


#def ine 


c_p_lCE 


2009.0 


/* 


J/ (kg K) 


specific heat capacity of ice */ 


#def ine 


rho_lCE 


910.0 


/* 


kg/(m-3) 


density of ice */ 


#def ine 


k_lCE 


2.10 


/* 


J/(m K s) 


= W/(m K) thermal conductivity of ice */ 


#def ine 


c_p_BRdef ault 


1000.0 


/* 


J/ (kg K) 


specific heat capacity of bedrock */ 


#def ine 


rho_BRdef ault 


3300.0 


/* 


kg/(m-3) 


density of bedrock */ 


#def ine 


k_BRdef ault 


3.0 


/* 


J/(m K s) 


= W/(m K) thermal conductivity of bedrock */ 


#def ine 


HO 


3000.0 


/* 


m */ 




#def ine 


BO 


1000.0 


/* 


m */ 




#def ine 


Ts 


223.15 


/* 


m */ 




#def ine 


G 


0.042 


/* 


W/(m-2) */ 


#def ine 


phi 


0.0125 


/* 


K/m */ 
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#define Nsum 30 /* number of terms in eigenf unction expansion; the exact 

solution is deliberately chosen to have finite expansion */ 



int exactK(const double t, const double z, double *TT, bool bedrocklslce) { 

int k; 

bool belowBO; 

double ZZ, P, alpha, lambda, beta, gamma, XkSQR, Xk, theta, Ck, 11, 12, aH, bB, ml, mR; 
double c_p_BR, rho.BR, k_BR; 

/* following constants were produced by calling print_alpha_k(30) (below) */ 
double alf[Nsum] = {3.350087528822397e-04, 1 . 114576827617396e-03 , 1 . 953590840303518e-03, 
2 . 684088585781064e-03 , 3 . 371114869333445e-03 , 4 . 189442265117592e-03 , 
5 . 008367405382524e-03 , 5 . 696044031764593e-03 , 6 . 425563506942886e-03 , 
7 . 264372872913219e-03 , 8 . 044853066396166e-03 , 8 . 714877612414516e-03 , 
9 .493529164160654e-03 , 1 . 033273985210279e-02 , 1 . 106421822502108e-02 , 
1 . 175060460132703e-02 , 1 . 256832682090360e-02 , 1 . 338784224692084e-02 , 
1 . 407617951778051e-02 , 1 . 480472324161026e-02 , 1 . 564331999062109e-02 , 
1 . 642470780103220e-02 , 1 . 709475346624607e-02 , 1 . 787248418996684e-02 , 
1 . 871188358061674e-02 , 1 . 944434477688470e-02 , 2 . 013010181370026e-02 , 
2 . 094721145334310e-02 , 2 . 176730968036079e-02 , 2 . 245631776169424e-02} ; 

if (bedrocklslce) { 
c_p_BR = c_p_lCE; 

rho_BR = rho_ICE; 
k_BR = k_ICE; 

for (k = 0; k < Nsum; k++) { /* overwrite alpha_k with ice-meets-ice values; see preprint */ 
alf [k] = (2.0 * k + 1.0) * pi / (2.0 * (HO + BO)); 

} 

} else { 

c_p_BR = c_p_BRdefault; 
rho.BR = rho.BRdef ault ; 
k_BR = k.BRdefault; 

} 

if (z > HO) { 
*TT = Ts; 
return 0; 

} 

belowBO = (z < -BO) ; 

ZZ = sqrt ( (rho.BR * c_p_BR * k_ICE) / (rho.lCE * c_p_ICE * k_BR) ) ; 
ml = (G / k_ICE) - phi; mR = (G / k_BR) - phi; 

/* DEBUG: printfC'ZZ = '/.lOe, ml = '/.lOe, mR = 7,10e\n" , ZZ,ml,mR); */ 
*TT =0.0; 

for (k = Nsum-l; k >= 0; k— ) { 

/* constants only having to do with eigenf unctions ; theta = theta_k(z) is the 

normalized eigenfunction */ 
alpha = alf [k] ; 
beta = ZZ * alpha; 

gamma = sin(alpha * HO) / cos (beta * BO); 

XkSQR = (rho.BR * c_p_BR * gamma * gamma * BO + rho.lCE * c_p_lCE * HO) / 2.0; 
Xk = sqrt (XkSQR) ; 

theta = ( (z >= 0) ? sin(alpha * (HO - z)) : gamma * cos(beta * (BO + z)) ) / Xk; 
lambda = (k_ICE * alpha * alpha) / (rho.ICE * c_p_lCE) ; 

/* DEBUG: printfC'k = "/.3d: alpha = "/.lOe, Xk = '/Me, theta = '/.lOe, lambda = •/.10e,\n", 

k, alpha, Xk, theta, lambda) ; */ 
/* constants involved in computing the expansion coefficients */ 
aH = alpha * HO; bB = beta * BO; 



EXACT TEMPERATURES IN ICE AND BEDROCK 



11 = - ml * (sin(aH) - aH * cos(aH)) / (alpha * alpha); 

12 = mR * (cos(bB) - 1.0 + bB * siii(bB)) / (beta * beta) 

- (BO * mR + HO * ml) * sin(bB) / beta; 
Ck = (rho.ICE * c_p_ICE * II + rho.BR * c_p_BR * gamma * 12) / Xk; 
/* add the term to the expansion */ 
*TT += Ck * exp(- lambda * t) * theta; 

/* DEBUG: printfC 11 = '/.lOe, 12 = '/.lOe, Ck = '/.lOe, term = "/.lOfSn", 

Il,I2,Ck, Ck * exp(- lambda * t) * theta ); */ 

} 

P = (z >= 0) ? (z / k.ICE) - (HO / k_ICE) : (z / k_BR) - (HO / k.ICE) ; 
*TT += Ts - G * P; 

return ( (belowBO) ? 1 : 0) ; 

} 



#define CQMPUTE.ALPHA 
#if CQMPUTE.ALPHA 

#define ALPHA.RELTOL l.Oe-14 
#define ITER.MAXED.OUT 999 

/* parameters needed for root problem: */ 
struct coscross_params { 

double Afrac, HZBsum, HZBdiff; 

}; 

/* the root problem is to make this function zero: */ 
double coscross (double alpha, void *params) { 

struct coscross_params *p = (struct coscross_params *) params; 

return cos (p->HZBsum * alpha) - p->Afrac * cos(p->HZBdif f * alpha); 

} 

/* compute the first N roots alpha_k of the equation 

((A-1)/(A+1)) cos((H - Z B) alpha) = cos((H + Z B) alpha) 
where H and B are heights and A, Z are defined In terms of material 
constants */ 

int print_alpha_k(const int N) { 

int status, iter, k, max_iter = 200; 
double Z, A; 

double alpha, alpha_lo, alpha_hi, temp_lo; 
const gsl_root_f solver_type *solvT; 
gsl_root_f solver *solv; 
gsl_f unction F; 

struct coscross_params params ; 

Z = sqrt((rho.BR * c_p_BR * k_ICE) / (rho_ICE * c_p_ICE * k_BR) ) ; 

A = (k_BR / k.ICE) * Z; 

params. Afrac = (A - 1 .0) / (A + 1 .0) ; 

params . HZBsum = HO + Z * BO; 

params . HZBdiff = HO - Z * BO; 

F. function = ftcoscross; 
F. params = ftparams; 

solvT = gsl_root_f solver _brent; // faster than bisection but still bracketing 
solv = gsl_root_f solver_alloc(solvT) ; 
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ED BUELER 



for (k = 0; k < N; k++) { 

// these numbers bracket exactly one solution 
alpha_lo = (double (k) * pi) / params.HZBsum; 
alpha_hi = (double (k + 1) * pi) / params.HZBsum; 
gsl_root_fsolver_set(solv, &F, alpha.lo, alpha.hi) ; 



iter = 0; 
do { 
iter++; 

status = gsl_root_f solver_iterate (solv) ; 

alpha = gsl_root_f solver _root (solv) ; 

alpha_lo = gsl_root_f solver_x_lower(solv) ; 

alpha_hi = gsl_root_f solver_x_upper(solv) ; 

temp.lo = (alpha.lo > 0) ? alpha.lo : (alpha_hi/2.0) ; 

status = gsl_root_test_interval (temp_lo , alpha_hi, 0, ALPHA_RELTOL) ; 
} while ((status == GSL.CONTINUE) &fe (iter < max.iter)); 
if (iter >= max_iter) { 

printf ("!!! ERROR : root finding iteration reached maximum iterations; QUITING!\n 

return ITER_MAXED_OUT; 

} 

printf ("y.ig. 15e,\n" .alpha) ; 

/* DEBUG: printf ("'/.le. 15e (in orig bracket ['/.IS. 15e ,7,19. 15e] )\n" .alpha, 

(double(k) * pi) / params.HZBsum, (double(k+l) * pi) / params.HZBsum); 

} 



gsl_root_f solver _free (solv) ; 
return status; 

} 

#endif /* COMPUTE.ALPHA */ 



